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Abstract 


Transformed  flux-form  semi-Lagrangian  (TFSL)  scheme  provides  stable  and  accurate 
algorithm  to  solve  the  advection-diffusion  equation.  Different  from  the  existing  flux-form  semi- 
Lagrangian  schemes,  the  flux  at  boundary  of  grid  cell  is  treated  as  the  temporal  mean  between 
present  and  next  time  steps.  After  the  temporal-spatial  transformation  using  the  characteristic-line 
concept,  the  temporal  integration  of  the  flux  from  present  to  next  time  step  becomes  the  spatial 
integration  of  the  flux  at  the  present  time  step.  This  scheme  is  always  stable  even  for  large  Courant 
numbers  (>  20)  with  the  second  order  accuracy  in  both  time  and  space.  For  the  Courant  number  not 
larger  than  0.5,  the  TFSL  scheme  reduces  to  the  Lax-Wendroff  scheme.  The  capability  of  the  TFSL 
scheme  is  demonstrated  by  the  simulation  of  the  equatorial  Rossby-soliton  propagation. 

Key  Words:  TFSL  scheme,  flux-form  semi-Lagrangian  scheme,  characteristic  line,  advection- 
diffusion  equation,  finite  volume,  conservative  finite  difference,  Courant  number,  equatorial 
Rossby  soliton. 
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1.  Introduction 


Numerical  approaches  in  atmospheric  and  oceanic  modeling  inevitably  introduce  diffusion 
(or  dissipation)  and  dispersion  into  the  approximate  solution.  From  a  physical  point  of  view, 
advection  of  a  passive  tracer  is  the  simple  transition  of  a  quantity.  Therefore,  dispersion,  the 
propagation  of  different  spatial  scales  at  different  phase  speed,  and  diffusion  are  processes  that  are 
aliens  to  the  process  that  is  being  modeled  (Chu  and  Fan  1998,  1999).  As  applied  to  constituent 
advection  problem,  these  numerical  artifacts  manifest  themselves  as  nonphysical  mixing  by 
numerical  diffusion,  nonphysical  highs  and  lows  in  the  constituent  field  caused  by  dispersion,  and 
nonphysical  tracer  spectra  caused  by  the  trapping  of  tracer  in  nonpropagating  small  spatial  scales 
(Rood  1987).  The  less  the  numerical  diffusion  and  dispersion  errors,  the  better  the  model 
performance  is. 

Propagation  of  a  Rossby  soliton  on  an  equatorial  beta-plane  is  treated  as  an  asymptotic 
solution,  which  exists  to  the  inviscid,  nonlinear  shallow  water  equations.  In  principle,  the  soliton 
propagates  to  the  west  at  fixed  phase  speed,  without  change  of  shape.  Since  the  uniform 
propagation  and  shape  preservation  of  the  soliton  are  achieved  through  a  delicate  balance  between 
linear  wave  dynamics  and  nonlinearity.  In  other  words,  the  Rossby  soliton  is  non-diffusive  and 
non-dispersive  (Boyd  1980),  which  makes  it  a  perfect  test  case  for  verification  of  numerical 
schemes  in  ocean  models  since  any  diffusion  and  dispersion  in  the  numerical  solution  of  the 
Rossby  soliton  are  computational  errors.  Interested  readers  are  referred  to  the  website: 
http://marine.rutgers.edu/po/index.php?model=test-problems. 

In  this  study,  we  first  show  instability  and  large  diffusion  and  dispersion  errors  in  numerical 
solution  of  the  Rossby  soliton  using  the  existing  schemes  such  as  the  upwind,  central,  and  Lax- 
Wendroff  schemes.  Then,  we  present  a  transformed  flux-form  semi-Lagrangian  (TFSL)  scheme. 
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which  has  explicit  form  and  much  less  diffusion  and  dispersion  errors.  The  numerieal  solution  of 
the  Rossby  soliton  exists  even  for  large  Courant  numbers. 

The  rest  of  paper  is  organized  as  follows.  Seetion  2  deseribes  the  equatorial  Rossby  soliton 
and  its  usefulness  for  the  oeean  model  verifieation.  Seetion  3  shows  the  failure  of  the  three  existing 
sehemes  (upwind,  central,  and  Lax-Wendroff)  in  simulating  the  equatorial  Rossby  soliton.  Seetion 
4  introduees  the  TFSL  seheme.  Seetion  5  derives  the  analytieal  form  of  the  amplifieation  faetor  of 
the  TFSL-seheme  and  shows  that  this  faetor  does  not  larger  than  1  for  large  Courant  number  sueh 
as  20.  Section  6  shows  the  eapability  of  the  TFSL-scheme  in  simulating  the  equatorial  Rossby 
soliton.  Seetion  7  presents  the  eonelusions. 

2.  Rossby  Soliton 


Let  Q  be  the  angular  frequeney  of  earth’s  rotation  and  R  be  the  earth  radius,  and  let  (x,  y)  be 
the  spatial  eoordinates  with  unit  veetors  (i,  j)  and  t  be  the  time.  Consider  a  single  layer  of 
homogeneous  oeean  layer  with  depth  of  H.  The  Lamb’s  parameter  is  defined  by 


E  = 


gH 


(1) 


where  g  is  the  gravitational  aeeeleration.  The  length  and  time  are  nondimensionalized  by 


L  = 


R 

771/4  ’ 


71/4 


T  = 


2Q 


(2) 


Let  (x,  y)  be  the  non-dimensional  Cartesian  eoordinates,  {u,  v)  be  the  non-dimensional  veloeity 
eomponents  in  the  meridional  and  latitudinal  directions,  and  (p  be  the  non-dimensional  surfaee 
elevation.  After  defining 

s  =  x-ct,  (3) 

and  transforming  the  nonlinear  shallow  water  wave  equations  into  a  frame  of  referenee  moving 
with  the  linear  wave,  the  flow  variables  {u,v,(/))  for  the  mode-1  ean  be  represented  by  (Boyd  1980) 
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u{s,y,t)  = 


(6/ -9) 


r]{s,t)  exp 


2 


(4a) 


v(.,v.0  =  2v^exp 

OS 


2 


(4b) 


<l>{s,y,t)  = 


(6/ +3) 


77(5,0  exp 


2 


(4c) 


and  the  variable  77(5,0  satisfies 


=  7i  =1-5366,  /,  =0.098765, 

ft  OS  OS 

which  is  the  Korteweg-de  Vries  (KDV)  equation  with  the  exact  solution, 

77(5,t)  =  ^sech"r5(5  +  /750)],  ^  =  0.77250  5  =  0.394,  /7  =  0.395, 


(5) 


(6) 


Substitution  of  the  exact  solution  (6)  into  the  third  term  in  the  lefthand  side  of  (5)  leads  to 

dt  ds 


(V) 


5(5,0  =  8yt5^/2{3sech'‘[5(5  +  /750)]  tanh[5(5  + 
-  sec  h^  [5(5  +  juB^t)]  tanh[5(5  +  juB^t)} } , 


(8) 


where  S  is  treated  as  a  source/sink  term.  Evidently  equation  (7)  has  the  analytical  solution  (6),  and 
therefore  it  can  be  used  to  verify  the  stability  of  the  numerical  schemes  since  the  diffusion  term  has 
been  changed  into  the  given  source/sink  term. 

3.  Several  Existing  Schemes 


To  solve  equation  (7)  numerically,  the  fluid  is  assumed  to  occupy  equatorial  region  around 
the  earth.  The  zonal  direction  is  discretized  into  120  cells  (i.e.,  resolution  at  3°  longitude).  The 
depth  of  fluid  is  set  up  as  100  m.  The  increment  A5  is  given  by 
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(9) 


A5  = 


IttR 

1201 


«  0.256. 


The  time  step  is  denoted  by  At .  Equation  (7)  ean  be  discretized  by  the  commonly  used  upwind 
scheme, 


Vi  =Vi  +C^{r|.^,-r|.)  +  S,^t , 


(10) 


the  central  scheme, 


n 


+  c;  + s';^t , 


and  the  Lax-Wendroff  scheme. 


vf =v:-^  ivL  -  vu ) + ^  (7;;.  -  2//; + vu  )+s:^t 


(11) 


(12) 


where  the  superscript  and  subscript  denote  the  time  step  and  the  horizontal  grid, 

V"=V{Xi,t„), 

and 

_  fxVl^t 

‘  Ay 


(13) 


(14) 


In  order  to  compare  the  difference  between  numerical  and  exact  solutions  (westward 
propagating  Rossby  soliton),  the  zonal  equatorial  strip  is  assumed  infinitely  long.  When  the  Rossby 
soliton  travels  over  nx  120  cells,  it  goes  around  the  earth  once  n  times  (called  n  cycles).  The  exact 
solution  at  t  =  0  is  taken  as  the  initial  condition, 

77(5,0)  =  sec /2^(&) ,  (15) 

with  5  =  0  denoting  0°  longitude. 
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The  three  differenee  equations  (10)-(12)  are  solved  numerieally  from  the  initial  eondition 
(15)  representing  the  upwind,  eentral,  and  Lax-Wendroff  sehemes  (Lax  and  Wendroff  1960)  with 
varying  At  at  eaeh  time  step  for  a  given  Courant  number  (C  =  0.75), 

A  CA5 

M  = - ^ — r- 

/imax(|;7;|) 

After  obtaining  the  numerieal  solution, //(x^t^),  substituting  it  into  (4o)  yields  (/){s^,y .  The 

aeeuraey  of  the  sehemes  ean  be  verified  through  their  eapability  in  predieting  the  westward 
propagation  of  the  Rossby  soliton.  To  do  so,  the  surfaee  elevation  is  plotted  with 

eontour  values  of  2.13,  4.26,  6.4,  8.53,  10.66,  12,79,  14.93,  and  17.06  em.  All  the  numerieal 
sehemes  greatly  distort  the  Rossby  soliton  (Fig.  1).  The  numerieal  solution  diverges  at  60°45’W 
using  the  upwind  seheme,  and  38‘’45’W  using  the  eentral  seheme.  The  numerieal  solution  does  not 
diverge  using  the  Lax-Wendroff  seheme,  however,  the  solution  is  totally  different  from  the 
analytieal  solution  after  propagating  one  eyele  around  the  earth  (eomparing  Fig.  Id  to  Fig.  la). 

4.  TFSL-Scheme 


4.1.  Semi-Lagrangian  Method 

Consider  the  adveetion  of  a  passive  sealar  by  the  veloeity  u(x,  t).  The  Eulerian 

formulation  of  this  is 


D(l) 

Dt 


d(j) 

dt 


+  u»V^  = 


(16) 


where  x  is  the  position  veetor,  DIDt  denotes  the  material  derivative,  while  the  Lagrangian 
eounterpart  is 


d(l>p  _ 
dt 


^  =  u(x^,0 


(17) 


7 


where  the  subseript  shows  the  fluid  partiele  in  Lagrangian  sense.  Although  (16)  and  (17)  earry 
the  same  physieal  information,  their  diseretization  and  numerieal  implementation  is  different:  (16) 
is  diseretized  on  an  Eulerian  grid  with  a  finite  number  of  grid  points  and  then  time-advaneed,  while 
(17)  is  integrated  for  a  finite  number  of  fluid  partieles. 

Semi-Lagrangian  methods  eombine  both  Eulerian  and  Eagrangian  points  of  view:  the  sealar 
field  is  diseretized  on  an  Eulerian  grid,  but  is  advaneed  in  time  using  (17).  The  key  element  in 
aeeomplishing  this  is  the  identifieation  of  eaeh  grid  point  x,  as  the  arrival  point,  for  instanee,  at 
t  +  At ,  of  a  partiele  originating  from  x*  at  time  t.  The  algorithm  has  three  steps:  (a)  The  partiele 

assoeiated  with  eaeh  grid  point  x,  at  time  t  +  At  is  traeed  baek  to  its  loeation  x*  at  time  t, 

^  /•  t+At 

X.  =x.-J  u(r)Jr;  (18) 

(b)  The  sealar  value  at  (x*,  t)  is  obtained  by  interpolating  the  known  values  at  neighboring  grid 
points, 

<l){x],t)  =  P[<l){{%j^},t],  (19) 

where  P  is  any  interpolation  operator  and  {  x  }  denotes  the  set  of  interpolation  points  assoeiated 
with  X* ,  for  example,  the  nodes  of  the  eell  eontaining  x* ;  (e)  Einally,  the  sealar  is  updated, 

(;^(x.,t  +  At)  =  ^zi(x*,t)  +  5'.At.  (20) 

Thus,  the  main  issues  of  the  semi-Eagrangian  method  are  the  baekward  integration  in  step  (a)  and 
the  interpolation  in  step  (b). 

4.2.  Flux  Form 

Equation  (16)  ean  be  rewritten  in  the  flux  form  with  inelusion  of  diffusion, 

|^  =  V.F  +  A,  F  =  -u^/i  +  xV^,  (21) 
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where  /eis  the  diffusion  eoeffieient.  Let  the  dependent  variable  be  defined  on  the  spaee  Q 

0<x<L^,  0<  y<  Ly,  0<  z<L^. 

with  {Lx,  Ly,  Lz)  the  lengths  in  (x,  y,  z)  direetions.  Let 

Ax  =  —,  Ay  =  —,  Az  =  — 

K  N, 

be  the  uniform  spatial  inerements  with  {Nx  +  Ny+  l,A^z+l)  the  grid  numbers.  Integrating  (21) 
for  the  finite  volume, 

^ijk  =  [^,-1/2  ^  ^  ^  yj-m  Lv+1/2.  Vi/2  ^  ^  ^  ^in/2] . 

_  Ax  _,Ay 

^i±l/2  —  —  2  ’  y  y±l/2  2  ’  ^k±l/l  —  ^k—  2  ’ 


from  to  +  At ,  we  obtain  the  finite  differenee  equation  of  the  flux-averaged  transport, 


^ij,k  ^i,j,k  _  ^  i+ll2,j,k  ^  i-l/2J,k  ,  ^  i,j+l/2,k  ^  ij-l/2,k  ,  Hij,k+\I2  Hij,k-\I2  ^  ^ 

+  +  +  '^ij,k  ’ 


At 


Ax 


Al 


Az 


(22) 


where  {F,  G,  H)  are  eomponents  of  the  veetor  F,  and 


-I  ^n+l 

f  Ydt, 

At  J 


(23) 


represents  the  temporal  average  (from  t„  to  t„+i).  The  tilde  represents  the  volume  average  over  , 


^ijk  ~ 


^ - III  (jkdxdydz  . 


AxA_yAz 


(24a) 


'ijk 


The  hat  represents  the  eombined  volume  (over  )  and  temporal  average  (from  t„  to  t„+i) 


1 


AtAxAyAz 


''n+\ 

I  III  Sdxdydzdt . 


(24b) 


^ijk 


For  the  finite  volume  AQ^^.^ ,  the  flux  at  x  =  x._j/2  and  t  =  t„  is  ealeulated  by 
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F” 


1 


I  ii2j,k  AyAz  ^yj-i/2 


c^k+\i2  rJ 
Jz/._i/o  Jv 


dx 


t — ty. 


dydz . 


(25) 


■i-l/2 


To  solve  equation  (22)  numerically,  we  need  to  compute  the  temporally  integrated  fluxes, 

«+l)  «+1)  «+l)  n+\)  «+l)  «+l) 

M/2J,k’  i-V2J,k’  /J+1/2,F  /J-l/2,i’  ^  i,J,k-l/2  ' 

If  these  fluxes  are  computed  using  the  semi-Lagrangian  method,  it  is  called  the  flux-form  semi- 
Lagrangian  scheme  (Lin  and  Rood  1996). 

4,3,  Transformation  of  Temporal  into  Spatial  Mean 

For  simplicity  and  no  loss  of  generality,  we  consider  one  dimensional  problem  of  (22) 
without  source/sink  term  (i.e.,  5,..^  =  0), 


7n+l  7n  77(«,n+l)  „{”."+0 

_Ku2  -Fi-V2 

At  Ax 


From  the  semi-Lagrangain  consideration,  we  have 


(26) 


ri7(«,«+l)  _  A/ 


Ax 


(27) 


Using  the  characteristic-line  concept,  the  flux  at  time  step  t„+i  and  location  x,.  1/2  can  be  transformed 
into  the  flux  at  time  step  t„  and  location  x,-.i/2-c  (Fig.  2), 

pn+\  pn  /28') 

1-1/2  1-1/2-C  ’ 

and  the  temporally  averaged  flux  F,-i/2  [similar  for  Fi+i/2  ]  can  be  transformed  into  spatial 
averaged  flux. 


(n,H+l)  1  rt„+ht  .  1  rXi_y2  N  , 


(29) 


where 
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(30) 


is  the  Courant  number.  Substitution  of  (30)  into  (29)  leads  to 


Z7«  I 


—(«,«+!) 
F i-\/ 2 


m-[ 

§  (i7«  +  p”  ] +  X  ( p”  P"  (p«  A.  p"  ) 

‘^1/2  p  1-1/2  ^  ^  i-1  )  i-k  ^  ^  i-k-l  ^mY  /-1/2-C  ) 


ifC<- 

2 


ifC>- 

2 


where 


1  1  1 
»>=  C--  .<^.-2=^.  <5.=^’ 


The  braeket  [  ]  represents  the  round-off  integer.  Similarly,  the  temporally  averaged  flux  at  the 
right  boundary  (x  =  x._^y2 ) 


pn  .  pn 

^i+\l2  ^i+\/2-C 


—(n,n+l) 

r  i+\i2 


rn—i 

^^1/2  i+1/2  ^  ^  i  Y  Z-I^k  i-k+\  ^  ^  i-k  ^mY  i+l-m  /+1/2-C  ) 


ifC<-^ 

2 


ifC>-^ 

2 


The  temporally  averaged  fluxes  and  (from  to  +  At )  is  transformed  into  the 

spatially  averaged  ones  over  multiple  grids  at  time  step  with  weights  of  d^i2,d^,...,5^.  If  the 
eharaeteristie  line  at  is  beyond  the  boundary,  the  boundary  eondition  ean  be  used  to  ealeulate 


FlY"  (Fig.  3), 


—(«,«+!) 

F2/2 


)(Fr+F-;2)+[c-0F, ■+/•■) 


where  is  the  boundary  value  between  F"  and  F^’ ,  and  is  interpolated  by 
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(35) 


K  = 


1 — 

I  2Cj 


‘  2C 


Substitution  of  (31)  and  (33)  into  the  difference  equation  (26)  leads  to 


in  F  \  jji  i 

<p.  =  ^^>.  + 


c<- 

2 


4 

T 


j(«,  +C,  -c, 


(36) 


c>i 

2 


which  is  called  the  Transformed  Flux-formed  Semi-Lagrangian  (TFSL)  scheme  for  the  advection- 
diffusion  equation  (21).  Here, 

D  =  C-m--. 

2 

For  C  <  1  /2,  the  TFSL  scheme  is  the  same  as  the  Lax-Wendroff  scheme.  Comparing  to  the 
central  difference  (CED),  the  TFSL-scheme  has  an  extra  positive  term, 


TSF-CED  =  ^(«(-,-2^;+C,), 


(37) 


for  C<l/2.  This  term  can  be  regarded  as  the  numerical  (positive)  diffusion  which  leads  to 
computational  stability.  Different  schemes  have  different  algorithms  to  compute  the  temporally 
averaged  fluxes  and  (from  to  C+At).  The  TFSL  scheme  has  second  order 

accuracy  in  both  time  and  space. 

5.  Stability  of  the  TFSL  Scheme 

Stability  of  numerical  schemes  is  an  important  issue  in  solving  the  advection  equation  (16). 
In  Section  3,  we  showed  the  instability  of  the  existing  schemes  (upwind,  central,  and  Lax- 
Wendroff).  To  determine  the  stability  of  the  TLSL  scheme  (36),  the  Lourier  series  expansion  is 
used.  Decay  or  growth  of  an  amplification  factor  indicates  whether  or  not  the  numerical  algorithm 
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is  stable  (von  Neumann  and  Riehtmyer  1950).  Assuming  that  at  any  time  step  4,  the  eompute 
solution  (f)"  is  the  sum  of  the  exact  solution  and  error  s" , 


(J)-  =  +  s" , 


(38) 


and  substituting  (38)  into  (36),  we  obtain 


i  i 


C<- 

2 


11  (  ^ 
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1 

—  ( s"  —  e"  'l h — I e"  —2e"+e'‘  'l 


(39) 


c>i 
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The  finite  mesh  function,  e"  ,  can  be  decomposed  into  a  Fourier  series, 

e:  =  I;  a;  exp(//^),  e  =  jnl  A,  (40) 

J=-N, 

with  /  =  ,  {a",  6)  being  the  amplitude  and  phase  angle  of  the  jth  harmonic.  Substituting  (40) 

into  (39)  yields 


=g{e,C)a'' 


(41) 


where 


g{0,c)  = 


l-C'(l-cos(9)-/Csin(9 

1  (  1  ^  f  ^  \ 

—  (l  -  cos  6)+  2  ~  ^  j  cos[(ot  - 1)  ^]  +  — cos  [(m  -  2) 

in  {mO)  + +  £)  -  £)^  j  sin  [(ot  - 1)  +  ^sin  [(m  -  2)  ^]| 


C<- 


(42) 


C>- 

2 


is  called  the  amplification  factor,  whose  magnitude  is  given  by 
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{o^c)\  = 


1  -  (1 -COS0)]  +C^  sin^  0 


ca 

2 
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(1  -COS0)^  + 


2  2 


2  A 


+  1  -  +  /)-Z)0  +  — 
2  j  4 


(43) 


+— (l  -COS0)^ 


cos 


(OT0)  +  |^a/)  jcos[(OT  -1)0]  +  ^cos[(ot  -2)^] 

cos (20) 
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^  1  D^'] 
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U  J  * 
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The  TFSL-scheme  is  computationally  stable  if  |g(^,C)|<l  and  computationally  unstable  if 

|g(^,C)|  >  1 .  Fig.  4  shows  that  |g(6',C)|  <  1  for  all  0  and  C  (larger  than  20),  which  implies  that 

the  TFSL-scheme  (36)  is  very  stable. 

6.  Simulating  the  Rossby  Soliton  Using  the  TFSL  Scheme 

The  TFSL-scheme  (36)  is  only  for  spatially  variant  and  temporally  invariant  u.  When  u  [or 
-fxV  in  (7)]  at  X/-1/2  varies  with  time  from  tn  to  tn+\,  concept  of  variant  characteristic  lines  can  be 

used  to  determine  m(x;.i/2,  t)  with  sub  time-steps  (<7^,2,  (between  tn  and  tn+\)  from  u{x,  U) 

at  grid  points  (x/.i,  . . .,  x,.m,  Xi*),  and  for  w  >  0  the  time  from  the  left  neighboring  grid  x ,.[^+1]  to  x^k 
is  given  by  (Fig.  5), 


_  ^'r‘  dx  _  Ax  r  dz 

^  J  ij(  ir  t  \  )/”  ,  1 


At  ln(l-C,_,)_  At  ,  C.  ,  Ch 


n(x,tj  Q 


C, 


-k 


n+lldL  +  lJdL  +  _\ 

2  3 


(44) 


k  =  0.5,  1, 
where 


n  n 


^i-k  ^i-[k+\] 


'^i-k 


_  ul,At 

Ax 


(45) 


The  parameter  C^k  is  the  Courant  number  for  sub  time  steps.  A  formula  similar  to  (44)  can  be 
obtained  for  m  <  0  (using  the  right  neighboring  grid).  The  temporally  averaged  fluxes  from  t„  to 
t„  +  At  can  be  calculated  by  (taking  F]*:"}" [see  (31)]  as  the  example) 


14 


p(n,n+\) 

i-Ml 


5t^ 


n 


{K 


m+l 


m-\ 


iLA+st 


(46) 


Equation  (7)  for  the  Rossby  soliton  is  discretized  using  the  flux  form, 


7,  -7, 

At 


p{n,n+\)  _  p{n,n+\) 
^i+l/2  ^i-Ml  I  ^ 

As 


(47) 


where  S-  is  the  temporally-spatially  averaged  source  term 


= 


1 


AtAs 


‘n+l 

i  i  S(s,t)dsdt . 


(48) 


with  S(s,  t)  given  by  (8).  The  difference  equation  (47)  is  solved  numerically  from  the  initial 
condition  (15)  using  the  TFSL-scheme.  To  compare  with  the  existing  schemes  (upwind,  central, 
and  Lax-Wendroff  schemes),  the  Courant  number  is  set  to  0.75.  After  the  numerical  solution 
ri{x.,t„)  is  obtained,  substituting  it  into  (4c)  yields  (/){s.,y j,tn)  ^  shown  in  Fig.  6.  Note  that  the 

three  existing  schemes  (upwind,  central,  and  Fax-Wendroff)  are  all  unstable  (Fig.  1),  but  the  TFSF- 
scheme  is  stable.  After  propagating  westward  around  the  earth  the  numerical  Rossby  soliton  (using 
the  TFSF  scheme)  shows  almost  non-difussive  and  non-dispersive. 

To  show  the  quality  of  the  TFSF-scheme,  the  difference  equation  (47)  is  integrated  for  C  = 
1.5  for  a  long  time  period  corresponding  to  the  Rossby  soliton  propagates  westward  around  the 
earth  5  times.  The  solution  (l){s^,y j,t^)  is  stable  all  the  time  (Fig.  7).  The  relative  root-mean-square 

error  (rrmse). 


rrmse(t)  = 


1  V,  A'v 


(nwm) 


{s.,y.,t)-(l)^"^\s„yj,t)^ 


max  (j)^‘‘^\s„y 


(49) 
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is  calculated  to  illustrate  the  aceuracy  of  the  TFSL  scheme.  Table  1  shows  RRMSE  at  the  end  of 
first  five  cycles  around  the  earth.  The  error  varies  from  2.66%  for  the  first  cycle  to  3.53%  for  the 
fifth  cycle. 

7.  Conclusions 

(1)  This  study  shows  that  TFSL  scheme  is  a  promising  stable  and  accurate  scheme  for 
solving  the  adveetion-diffusion  equation.  Magnitude  of  the  amplification  factor  does  not  exceed  1 
for  large  Courant  number  (e.g.,  for  C  =  20).  The  Fourier  analysis  shows  that  the  TFSL  seheme  has 
second-order  aecuracy  in  time  and  space.  Computational  stability  and  higher  aceuraey  than  the 
widely  used  schemes  (eentral,  upwind,  and  Lax-Wendroff)  makes  this  scheme  useful  in  oeean 
modeling,  computational  fluid  dynamics,  and  numerical  weather  prediction. 

(2)  Several  major  features  distinguish  the  TLSL  seheme  from  the  existing  schemes,  both 
Lulerian  and  semi-Lagrangian.  Lirst,  the  flux  (F)  at  the  boundary  of  eaeh  grid  eell  is  eomputed  not 
from  a  single  time  step  (present  or  next)  but  from  temporal  integration  from  present  to  next  time 
step.  Second,  this  temporal  integration  is  transformed  into  spatial  integration  at  the  present  time 
step  using  the  characteristic  line  method. 

(3)  The  equatorial  Rossby  soliton  is  used  to  test  the  eapability  of  the  TLSL  scheme  since  it 
has  exact  solution.  The  equation  is  solved  numerieally  from  the  soliton  initially  loeated  at  the 
equator  and  0“  longitude  with  the  overall  Courant  number  of  0.75.  The  existing  numerical 
schemes  greatly  distort  the  Rossby  soliton  and  diverge  as  it  propagates:  60°45’W  (upwind), 
38“45’W  (eentral),  and  totally  distorted  after  one  cyele  around  the  earth  (Lax-Wendroff).  However, 
the  TLSL  seheme  does  not  distort  the  Rossby  soliton  and  converge  as  it  propagates  many  cyeles 
around  the  earth. 
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(4)  Future  studies  include  applying  TFSL  scheme  to  non-uniform  grid  systems  as  well  as 


designing  higher  order  TFSL  schemes. 
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Table  1.  RRMSE  of  the  surfaee  elevation  predieted  using  the  TFSL-scheme  after  the  first  five 
cycles  around  the  earth. 


Cycle 

1 

2 

3 

4 

5 

RRMSE  (%) 

2.66 

2.86 

3.00 

3.22 

3.53 
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(a):  Analytic  Solution  with  Courant  Number=  0.75 


(d):  Lax-Wendroff  Scheme  One  Cycle 
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Fig.  1.  Surface  elevation  (l){s,y,t)oi  the  Rossby  solitons  obtained  from  (a)  exact  solution,  and 
numerical  integration  with  C  =  0.75  using  (b)  upwind  scheme,  (c)  central  scheme,  and  (d)  Lax- 
Wendroff  scheme.  Note  that  the  numerical  solution  diverges  at  60“45’W  using  the  upwind  scheme 
and  38°45’W  using  the  central  scheme.  The  numerical  solution  does  not  diverge  using  the  Lax- 
Wendroff  scheme,  however,  the  solution  is  totally  distorted  from  the  analytical  solution  after 
propagating  one  cycle  around  the  earth. 
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Fig.  2.  Temporally  varying  flux  at  the  boundary  Xi.1/2  from  to  +  At  is  transformed  into  spatially 
varying  flux  at  t  n  from  Xi.1/2  -  CAt  to  Xi.1/2  using  the  characteristic-line  concept. 
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CFL  number 


(a):  Analytic  Solution  with  Courant  Number=  0.75 


Fig.  6.  Surface  elevation  (l){s,y,t)oi  the  Rossby  soliton  obtained  from  (a)  exact  solution,  and 
numerical  integration  with  C  =  0.75  using  (b)  TFSL-scheme.  The  solutions  (l){s,y,t)  are  plotted  at 
four  time  instances  for  the  Rossby  soliton  (exact  solution)  westward  propagating  90°,  180  °,  270°, 
and  360°  (return  to  the  initial  location). 
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